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The dynamic structure factor S(k, ui) and the two-particle 
distribution function g(r, t) of ions in a Coulomb crystal are 
obtained in a closed analytic form using the harmonic lattice 
(HL) approximation which takes into account all processes 
of multi-phonon excitation and absorption. The static ra- 
dial two-particle distribution function g(r) is calculated for 
classical (T > Tiuj p , where uj p is the ion plasma frequency) 
and quantum (T <C hu> p ) body-centered cubic (bcc) crys- 
tals. The results for the classical crystal are in a very good 
agreement with extensive Monte Carlo (MC) calculations at 
1.5 < r/a < 7, where a is the ion-sphere radius. The HL 
Coulomb energy is calculated for classical and quantum bcc 
and face-centered cubic crystals, and anharmonic corrections 
are discussed. The inelastic part of the HL static structure 
factor S"(k), averaged over orientations of wave- vector k, is 
shown to contain pronounced singularities at Bragg diffrac- 
tion positions. The type of the singularities is different in 
classical and quantum cases. The HL method can serve as a 
useful tool complementary to MC and other numerical meth- 
ods. 

PACS numbers: 52.25.Zb 



I. INTRODUCTION 

A model of a Coulomb crystal of point charges in 
a uniform neutralizing background of charges of oppo- 
site sign is widely used in various branches of physics. 
The model was originally proposed by Wigner Q who 
showed that zero-temperature electron gas immersed into 
uniform background of positive charges crystallizes into 
body-centered cubic (bcc) Coulomb crystal at sufficiently 
low density. Since then the model has been used in solid 
state physics for describing electron- hole plasma (e.g., 
Ref . Q ) and in plasma physics for describing dusty plas- 
mas and ion plasmas in Penning traps (e.g., Ref. Q). 
Finally, Coulomb crystals of ions on almost uniform back- 
ground of degenerate electron gas are known to be formed 
in the cores of white dwarfs and the envelopes of neutron 
stars. Consequently, properties of Coulomb crystals are 
important for studying structure and evolution of these 
astrophysical objects (e.g., Ref. [Q). 

As classical examples of strongly coupled systems, the 



Coulomb crystals have been the subject of extensive stud- 
ies by various numerical methods, mostly by Monte Carlo 
(MC; e.g., ||, and references therein), and also by molec- 
ular dynamics (MD; e.g., Ref. Q), and path-integral 
Monte Carlo (PIMC; e.g, Ref. @). Although the results 
of these studies are very impressive, the numerical meth- 
ods are time consuming and require the most powerful 
computers. 

The aim of the present article is to draw attention to a 
simple analytic model of Coulomb crystals. It has been 
employed recently in Ref. Q in connection with trans- 
port properties of degenerate electrons in strongly cou- 
pled plasmas of ions. We will show that this model is a 
useful tool for studying static and dynamic properties of 
Coulomb crystals themselves. 



II. STRUCTURE FACTORS IN HARMONIC 
LATTICE APPROXIMATION 

For certainty, consider a Coulomb crystal of ions im- 
mersed in a uniform electron background. Let p{jc,t) = 
J2i — rj(t)) be the Heisenberg representation opera- 
tor of the ion number density, where f i (t) is the operator 
of the ith ion position. The spatial Fourier harmonics of 
the number density operator is Pk(t) = ^ e -tkT 1 (t) > ^ nc 

dynamic structure factor S(k, oj) of the charge density is 
defined as 
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where N is the number of ions in the system, n is the 
ion number density, (. . .)t means canonical averaging at 
temperature T, and the last term takes into account con- 
tribution from the neutralizing background. 

The above definition is equally valid for liquid and solid 
states of the ion system. In the solid regime, it is natural 
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to set Ti(t) = Rj + u-t(i), where is a lattice vector, 
and Uj(t) is an operator of ion displacement from R^. 
Accordingly, 



where n v — (e Zv — 1) 1 is the mean number of phonons 
in a mode v. The brackets 
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The main subject of the present paper is to discuss 
the harmonic lattice (HL) model which consists in re- 
placing the canonical averaging, (. . .}t, based on the ex- 
act Hamiltonian, by the averaging based on the corre- 
sponding oscillatory Hamiltonian which will be denoted 
as (. . .)to- In order to perform the latter averaging we 
expand Ui(i) in terms of phonon normal coordinates: 
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where m is the ion mass, v = (q, s), s = 1,2,3 enumer- 
ates phonon branches; q, e v , lo v are, respectively, phonon 
wavevector (in the first Brillouin zone) , polarization vec- 
tor, and frequency; b v and b\ refer to phonon annihilation 
and creation operators. The averaging over the oscilla- 
tory Hamiltonian, H — J2 u \h,uj v {b v b\ + i>l,b v ), reads 



(5) 



where n u is the number of phonons in a mode v, f{n v ) — 
e~ n " z " (1 — e~ z ") is the phonon density matrix in thermo- 
dynamic equilibrium, z„ = Hlu^/T, F nvnv is a diagonal 
matrix element of the operator F . Inserting Eq. (gj) into 
(||) we can perform the averaging (0) using the technique 
described, for instance, in Kittel |gf. 

The resulting structure factor 5(k, t) takes into ac- 
count absorption and emission of any number of phonons; 
it can be decomposed into the time-independent elastic 
(Bragg) part and the inelastic part, S(k, t) = <S"(k) + 
S"(k,t). The elastic part is f§: 



S'(k) = e- 2W ^ {2nfn 



G 



(6) 



where G is a reciprocal lattice vector; prime over the sum 
means that the G = term is excluded (that is done due 
to the presence of uniform electron background) . 

In Eq. (0) we have introduced the Debye- Waller factor, 
e -W(k) = Yexp(«k- u)), A 
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denote averaging over the phonon spectrum, which can 
be performed numerically, e.g., Ref. |10| . The integral 
on the rhs is meant to be taken over the first Brillouin 
zone. The latter equality in Eq. (Q) is exact at least 
for cubic crystals discussed below. For these crystals, 
W(k) = r^k 2 /6, where r% = (u 2 )to is the mean-squared 
ion displacement (e.g., P,Eof). 

The inelastic part of 5(k, t) (e.g., can be rewritten 

as 
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Eqs. ([ 
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and (H) result in the HL dynamical structure 
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where t is real and r = t — ih/(2T). 

Along with the HL model we will also use the simpli- 
fied model introduced in Ref. ||. It will be called HL1 
and its results will be labelled by the subscript '1'. It 
consists in replacing S"(k,t) given by Eq. (||) by a sim- 
plified expression S"(k,t) equal to the first term of the 
sum, R = 0: 



Si(k,t)=S'(k)+S'{(k,t), 
S'l(k,t)=e- 2W ^ (VW fc2 - 
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where v is defined by the equation v a p(0,t) — v{t)5 a p, 
which is the exact tensor structure for cubic crystals (see 
above). The accuracy of this approximation, as discussed 
in Ref. || , is good for evaluating the quantities obtained 
by integration over k (e.g., transport properties of degen- 
erate electrons in Coulomb crystals of ions). 



III. STATIC CASE. HL VERSUS MC 

In this section we compare our analytic models with 
MC simulations of Coulomb crystals. For this purpose 
we introduce the function 
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which may be called the static two particle radial distri- 
bution function. This function is the result of an angular 
and a translation average of the static two particle dis- 
tribution function. In this expression df2 r is the solid 
angle element in the direction of r. One can see that 
4irr 2 ng(r)dr is the ensemble averaged number of ions in 
a spherical shell of radius r and width dr centered at a 
given ion. Thus g(r) is just the quantity determined from 
MC simulations 0. 

First let us use the HL1 model. From Eqs. (j|) and JT^ ) 
we easily obtain gi(r) = g'ir) 4- g'{(r), where 
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Calculation of g" (r) in the HL model is more cumber- 
some. After integration over k = |k| and £l r the result 
can be written as 



9{r) = gx{r) 
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= 3(r + GRf/{Ar 2 T ), 
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where 7 = (r + aRfj,)/x, rj 

cosd, $ is an angle between k and R, x 2 = 4[r^/3 
(k a kpVap(R,0)/k 2 )], and df^k is the solid angle element 
in the direction of k. Therefore, we need to evaluate 
a rapidly converging lattice sum ( |l5| ) of 2D integrals in 
which x is known once the matrix elements v Q ^(R, 0) 
are calculated from Eq. (10). We have performed the 
integration over the first Brillouin zone required in Eq. 
( |i"o| ) using the 3D Gauss integration scheme described in 
Ref. @. 

The function g(r) depends on the lattice type and on 
two parameters: the classical ion coupling parameter T = 
Z 2 e 2 / (aT) and the quantum parameter 8 — Huj p /T that 
measures the importance of zero-point lattice vibrations. 
In this case Ze is the ion charge, a — (47rn/3)~ 1 / 3 is the 
ion sphere radius, and uj p = Ze \J Aim/m the ion plasma 
frequency. 

First consider a classical Coulomb crystal, 9^0, for 
which n v sa T / (htu u ). The functions g(r) calculated using 
the HL and HL1 models for body-centered cubic 
crystals at T = 180 and 800 are presented 
and |^. The pronounced peak structure corresponds to 
the bcc lattice vectors. These results are compared with 
extensive MC simulations. The MC method is described, 
e.g., in Ref. ||. The simulations have been done with 686 
particles over nearly 10 8 MC configurations. 
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FIG. 1. g(r) for a bcc Coulomb crystal at V = 180. 

One can observe a very good agreement of HL and 
MC results for both values of T at 1.5 < r/a < 7. The 
MC results for g(r) are limited to half the size of the 
basic cell containing the N charges due to the bias from 
particles in the image cells adjacent to the basic cell. 
For N = 686 the basic cell length is 14.2 a. Hence the 
MC g(r) results for this simulation are valid only out 
to r ps 7 a while g(r), given by the HL model, remains 



accurate as r — > 00. At small particle separations, 



< 

1.5 a, where g(r) becomes small, the HL g(r) deviates 
from the MC g(r). It is clear that the HL model cannot 
be reliable at these r, where strong Coulomb repulsion 
of two particles dominates, and the MC data (available 
down to r ^ 1.1 a) are more accurate. The HL1 model is 
quite satisfactory at r > 2.5 a, beyond the closest lattice 
peak. The HL model improves significantly HL1 at lower 
r. It is interesting that for T = 180 the HL1 model agrees 
slightly better with MC for the range 2.5 < r/a < 6 than 
the HL model does. With increasing T, however, the HL 
model comes into better agreement with MC at these r, 
although the difference between the HL and HL1 models 
becomes very small. This good agreement of the HL 
models with the MC simulations after the first peak of 
g(r) indicates that we have a very good description of 
Coulomb crystals for which the HL model may be used 
in place of MC simulations. 

The HL model enables one to analyse quantum effects. 
Figs, [l] and |^ exhibit also g(r) in the quantum regime 
at 8 = 10. Zero-point lattice vibrations tend to reduce 
lattice peaks. The simplicity of the implementation of 
the HL model in the quantum regime is remarkable given 
the complexity of direct numerical studies of the quantum 
effects by MC, PIMC or MD simulations (see, e.g., Ref. 

0). 
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FIG. 2. Same as in Fig. fU but at V = 800. 



IV. COULOMB ENERGY 

To get a deeper insight into the HL and HL1 models 
let us use them to calculate the electrostatic energy U of 
the crystal. Writing this energy as the sum of Coulomb 
energies of different pairs of ions complemented by the 
interaction energy of ions with the electron background 
and the Coulomb energy of the background itself, we ar- 
rive at the standard expression 
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where g(r) is given by Eq. (fl5|). Therefore, we can use 
the function g(r) calculated in Sect. 3 to analyse U. 
For the HL1 model from Eqs. (|T|) we get 
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where £ is the electrostatic Madelung constant [= 
-0.895929 for bcc, and -0.895873 for face-centered cu- 
bic (fee) lattice], and erfc(cc) is the complementary error 
function. The second line of this equation is obtained us- 
ing the formula for the Madelung constant derived with 
the Ewald method (see, e.g., Ref. Jl^| ) 
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where A is an arbitrary number. In the particular case 
of Eq. @ A = V3a/(2r T )- 



For the HL model, using Eq. (15), we have 
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consider the classical crystal at zero tempera- 
— > 0. Then rx — > 0, x — > 0, and we reproduce 



R 

First, 
ture, T 

the Madelung energy, U/N -> Ui/N C^ 2 e 2 /a. In the 
limit of small T both Ui/N and f//7V contain the main 
term that can be expanded in powers of T plus an ex- 
ponentially small term (non- analytic at T = 0). For the 
classical crystal at any T we have r^/a 2 = «_2/r, where 
U s = ((u>i//uJp) B )ph denotes a phonon spectrum moment 
(u_ 2 =12.973 for bcc and 12.143 for fee). 

The sum over R ^ in the last expression for U\ 
in Eq. ( |l7j ) is exponentially small. Thus the analytic 
part of U\ in the HL1 model is given only by two terms, 
Ui/(NT) = Cr + u_ 2 /2. We see that the HL1 model fails 
to reproduce correctly the harmonic part of the potential 
energy: u-2/2 appears instead of conventional 3/2. 

On the contrary, the expansion of UJ (NT) in the HL 
model, Eq. (|l9|), contains all powers of T. To analyse 
this expansion, let us take any term of the sum over R, 
and introduce a local coordinate frame with z-axis along 
R. Then 
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where <j> is an azimuthal angle of k in the adopted frame. 
Since x — > as T — > in the denominator of the exponent 
under the integral in Eq. (|l^), only a narrow interval of 
fj, in the vicinity of fj, — contributes, and we can extend 
the integration over /i to the interval from — 00 to +00. 
Furthermore, using the definition of x, Eq. (|l5|), we can 
rewrite x as 
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where v a p — v a p (R, 0) . Accordingly, we can treat e as 
small parameter and expand any integrand in Eq. jl9| ) 
in powers of e and further in powers of fi. This generates 
the expansion in powers of T. 

We have been able to evaluate three first terms of this 
expansion. In particular, the term linear in T contains 
the expression 
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where 2? a( 3 is the dynamical matrix. Combining this 
expression with r 2 ^/(2a 2 ) and taking into account that 
T^ape-uae-vf} = ^l/ujp (according to the basic equa- 
tion for the phonon spectrum) we see that the HL ex- 
pansion of the analytic part of U in powers of T is 
U/(NT) = Cr + 3/2 + SU T /(NT); it reproduces not 
only the Madelung term, but also the correct oscilla- 
tory term 3/2, and contains a higher-order contribution 
SU T /(NT) = Af L /T + A% L /T 2 + ... that can be called 
"anharmonic" contribution in the HL model. After some 
transformations the coefficients and are re- 

duced to the sums over R containing, respectively, bi- 
linear and triple products of v a $ (with integration over 
fi and (j) done analytically). Numerically the sums yield 
Af L = 10.64 and Af L = -62.4. 

The anharmonic terms occur since U, as given by Eq. 
( |l6| ) , includes exact Coulomb energy (without expanding 
the Coulomb potential in powers of ion displacements 
u). However, we use g(r) in the HL approximation and 
thus neglect the anharmonic contribution in ion-ion cor- 
relations. Therefore, the HL model does not include all 
anharmonic effects. 

Let us compare the HL calculation of SUt with the 
exact calculation of the first anharmonic term in the 
Coulomb energy of classical Coulomb crystals by Dubin 
@. The author studied the expansion SU T xact /{NT) = 

and expressed the first term 
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where U n /n\ is the nth term of the Taylor expansion of 
the Coulomb energy over ion displacements, while angu- 
lar brackets denote averaging with the harmonic Hamil- 
tonian H . According to Dubin Af xact =10.84 and 12.34 
for bec and fee crystals, respectively. (The same quantity 
was computed earlier by Nagara et al. |l4| ] who reported 
Af xact =10.9 for bcc.) 

It turns out that our SUt sums up a part of the 
infinite series of anharmonic corrections to the energy, 
denoted by Dubin as J^La^n) / '("Oj so that A^ h = 
r(f7 4 )/(24A/T), A$ L = r 2 ([/ 6 )/(6WT), etc. (The fact 
that this summation can be performed in a closed ana- 
lytic form was known from works on the so called self- 
consistent phonon approximation, e.g., |pl| and refer- 
ences therein.) Our numerical value for the bcc lattice 
Af h = 10.64 is very close to the value of T{U A )/{2ANT) 
reported by Dubin as w 10.69 (his Table 3) which 
confirms accuracy of both calculations. The fact that 
Af L = 10.64 is close to Af act = 10.84 for bcc is acci- 
dental (Dubin found r(£/ 3 2 )/(72ATT 2 ) « 21.53 for bcc). 



For instance, from the results of Ref. |l3j for fee one in- 
fers, A^ L ps 5.63 which differs strongly from the exact 
anharmonic coefficient ^4° xact = 12.34. 

Now let us set T = and analyse the quantum effects. 
We can expand Eqs. ( p"7j ) and ( |l9| ) in powers of Tt/cl. For 
T = the quantity tt tends to the rms amplitude of zero- 
point vibrations, Tt = \JZTiu_il (2muj p ), where u_i is 
another phonon spectrum moment (=2.7986 and 2.7198 
for bcc and fee, respectively). The expansion of U\/N 
gives CZ 2 e 2 /a + u_ihojp/4 plus small non-analytic terms. 



In the same manner as in Eq. fl22p we find that U/N 
C,Z 2 e 2 1 a + Suihujp/A + SUq/N. The second term gives 
half of the total (kinetic + potential) zero-point harmonic 
energy of a crystal, as required by the virial theorem for 
harmonic oscillator (ui =0.51139 and 0.51319 for bcc and 
fee, respectively), while the third term, SUq, represents 
zero-point anharmonic energy in the HL approximation. 

To make the above algebra less abstract let us esti- 
mate the accuracy of the HL model and the relative im- 
portance of the anharmonicity and quantum effects. In 
the classical case, taking T = 170 (close to the melt- 
ing value r m = 172 for bcc), we estimate the anhar- 
monic contribution to the total electrostatic energy as 
\SU T /U\ w ^f act /(|C|r 2 ) w 4.2 x 10" 4 and 4.8 x 10~ 4 
for bcc and fee, respectively. 

The relative error into U introduced by using the HL 
model is A!f act /(|C|r 3 ) » 5.7 x 10~ 5 for bcc (if we adopt 
an estimate of A!f act m 247 from the MD data on the full 
electrostatic energy presented in Table 5 of Ref. || ) and 
j^oxact _ A HLy^Q T 2j _ 2 .6 x 10~ 4 for fee. We see that 
Coulomb crystals can be regarded as highly harmonic, 
and the accuracy of the HL model is sufficient for many 
practical applications. Obviously, the accuracy becomes 
even better with decreasing T. The quantum effects can 
be more important (than the anharmonicity) in real sit- 
uations. Let us take 12 C matter at density p = 10 6 g 
cm~ 3 typical for the white dwarf cores or neutron star 
crusts. The quantum contribution into energy is mea- 
sured by the ratio 3uihui p / '(4|£|Z 2 e 2 / 'a) which is equal 
to 4.7 x 10~ 3 at given p (and grows with density as p 1 ^ 6 ). 

For completeness we mention that the compressibility 
of the electron background also contributes to the electro- 
static energy. The relative contribution in the degenerate 
electron case for 12 C at p = 10 6 g cm~ 3 is ~ 10 -2 (e.g., 
Ref. fll6| ). Another point is that the HL model takes into 
account zero-point lattice vibrations but neglects ion ex- 
change which becomes important at very high densities 
(e.g., Ref. |). 



V. STRUCTURE FACTORS 

Finally, it is tempting to use the HL model for an- 
alyzing the ion structure factors themselves. Con- 
sider the angle-averaged static structure factor S(k) = 
fdn k S(k,t = 0)/(4tt). For the Bragg part, from Eq. (§) 
we obtain the expression 
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containing delta-function singularities at k = G, lengths 
of reciprocal lattice vectors G. Direct HL calculation of 
S"(k) from Eq. (|9|) is complicated by the slow conver- 
gence of the sum and complex dependence of v a p on R. 
However, the main features of S"(k) can be understood 
from two approximations. First, in the HL1 model we 
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have v a/3 (0,0)k a kp = 2W{k), and S£(k) 
as shown by the dashed line in Fig. [|. 

The second, more realistic approximation will be 
called HL2 (and labelled by the subscript '2'). It 
consists in adopting a simplified tensor decomposition 
of v a p(R, 0) of the form u Oj g(R,0) = F(R)S a p + 
R a RpJ(R)/R 2 . If so, we can immediately take the fol- 
lowing integrals / dft R v aa (R, 0)/(4tt) = 3F(R) + J(R) 
and J dn R v al3 (R,O)R a R /(4wR 2 ) = F(R) + J(R) (as- 
suming summation over repeating tensor indices a and 
0). On the other hand, we can calculate the same inte- 
grals taking u aj a(R, 0) from Eq. ( pX| ) at t = 0. In this 
way we come to two linear equations for F(R) and J(R). 
Solving them, we obtain 
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where y = qR, and jo(y) and ji(y) are the spherical 
Bessel functions. Note that F(0)fc 2 = 2W(k), J(0) = 0. 

In the limit of large R the functions jo(qR) and ji(qR) 
in Eqs. ( |25| ) strongly oscillate which means that the main 
contribution into the phonon averaging (integration over 
q) comes from a small vicinity near the center of the Bril- 
louin zone. Among three branches of phonon vibrations 
in simple Coulomb crystals, two (s=l, 2) behave as trans- 
verse acoustic modes, while the third (s=3) behaves as 
a longitudinal optical mode (u> ~ u) p ) near the center of 
the Brillouin zone. Owing to the presence of ui^ 1 in the 
denominator of Eqs. (25), the main contribution at large 
R comes evidently from the acoustic modes. Thus we 
can neglect optical phonons and set u> — c s q for acoustic 
modes, where c s is the mean ion sound velocity. In the 
high-temperature classical limit, {n v + g) — ► T/(Kc s q). 
Then from Eqs. (|25|) at R — ► oo we approximately obtain 
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Our analysis shows that an appropriate value of cj~ 
for bcc lattice would be 67.85/(aw p ) 2 . From Eq. ( p6[ ) we 
see that F(R) and J(-R) decrease as R^ 1 with increasing 
i?. In the quantum limit > 1 we have (n v + 5) - * §; 
applying the same arguments we deduce that F, J ex. R~ 2 
as i? — > 00. 

Using Eq. (0) we have 
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(27) 



A number of the first terms of the sum, say for |R| < 
Rq, where Ro/a is sufficiently large, can be calculated 
exactly. To analyse the convergence of the sum over R 
at large R let us expand the exponential in the square 
brackets on the rhs. All the terms of the expansion which 
behave as R~ n with n > 2 lead to nicely convergent 
contributions to S'^k). The only problem is posed by 
the linear expansion term in the classical case. The tail 
of the sum, X)|r|>_r ' ^ or ^ s term can be regularized and 
calculated by the Ewald method (e.g., Ref. |fL2f ) with the 
following result 
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where Ei(— x) is the exponential integral, and A is a num- 
ber to be chosen in such a way the convergence of both 
infinite sums (over direct and reciprocal lattice vectors) 
be equally rapid. Letting A — * 00 we obtain a much more 
transparent, although slower convergent formula 
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This expression explicitly reveals logarithmic singular- 
ities at A; = G. They come from inelastic processes of 
one-phonon emission or absorption in the cases in which 
given wave vector k is close to a reciprocal lattice vec- 
tor G. To prove this statement let us perform Taylor 
expansions of both exponentials in angular brackets in 
Eq. (§). The one-phonon processes correspond to those 
expansion terms which contain products of one creation 
and one annihilation operator. Thus, in the one-phonon 
approximation S"(k, t = 0) reads 
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where the last summation is over phonon polarizations, 
q = k — G is the phonon wave vector which is the given 
wave vector k reduced into the first Brillouin zone by 
subtracting an appropriate reciprocal lattice vector G. 
In addition, in Eq. ( pp| ) we have introduced an over- 
all factor e - 2W ( k ) which comes from renormalization of 
the one-phonon probability associated with emission and 
absorption of any number of virtual phonons (e.g., Ref. 
||). Now let us assume that \k — G\a -C 1 and average 
Eq. (30) over orientations of k [integrate over dfik/(47r)]. 
One can easily see that the important contribution into 
the integral comes from a narrow cone tto aligned along 
G. Let So < 1 be the cone angle chosen is such a way 
that Gdoa <C 1, but G9q 3> \G — k\. Integrating within 
this cone, we can again adopt approximation of acoustic 
and longitudinal phonons and neglect the contribution 
of the latters. For simplicity, we also assume that the 
sound velocities of both acoustic branches are the same: 
w u = c s |k — G|. Then, in the classical limit we come to 
the integral of the type 
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which contains exactly the same logarithmic divergency 
we got in Eq. (^9|). Note that in the quantum limit we 
would have similar integral but with u> instead of to 2 in 
the denominator of the integrand. The integration would 
yield the expression proportional to \k — G\, i.e., the log- 
arithmic singularity would be replaced by a weaker kink- 
like feature. Therefore, the k = G features of the inelastic 
structure factor S"(k) in the quantum limit are expected 
to be less pronounced than in the classical limit but could 
be, nevertheless, quite visible. Actually, at any finite 
temperature, even deep in the quantum regime T <C hto p 



there are still phonons excited thermally near the very 
center of the Brillouin zone, where the energy of acous- 
tic phonons is smaller than temperature. Due to these 
phonons the logarithmic singularity always exists on top 
of the kink-like feature at T ^ 0. 

After this simplified consideration let us return to qual- 
itative analysis. We have calculated S^' (k) in the classical 
limit using the HL2 approximation as prescribed above 
and verified that the result is indeed independent of R 
(in the range from ~ 30a to 100a) and A. The resulting 
S'2 (k) is plotted in Fig. ||by the solid line. 
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FIG. 3. Inelastic part of the structure factor at T — 180 for 
classical bcc crystal. 

Thus, in a crystal, the inelastic part of the structure 
factor, S"(k), appears to be singular in addition to the 
Bragg (elastic) part S'(k). The singularities of S"(k) 
are weaker than the Bragg diffraction delta functions in 
S'(k); the positions of singularities of both types coin- 
cide. The pronounced shapes of the S"(k) peaks may, 
in principle, enable one to observe them experimentally. 
The structure factor S(k) in the Coulomb liquid (see, e.g., 
Ref. [0 and references therein) also contains significant 
but finite and regular humps associated with short-range 
order. This structure has been studied in detail by MC 
and other numerical methods. In contrast, the studies 
of singular structure factors in a crystal by MC or MD 
methods would be very complicated. Luckily, they can 
be explored by the HL model. 

Finally, it is instructive to compare the behavior of 
S"(k) at small k in the HL1 and HL2 models. It is easy 
to see that the main contribution to inelastic scatter- 
ing at these k comes from one-phonon normal processes 
[with q=k in Eq. ©]. At these k the HL2 S'^ik) co- 
incides with the one-phonon S" b (k) and with the static 
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structure factor of Coulomb liquid (at the same T) and 
reproduces correct hydrodynamic limit S(k) cx k 2 . 
The HL1 model, on the contrary, overestimates the im- 
portance of the normal processes. 

Let us mention that we have also used the HL2 model 
to calculate g(r). HL2 appears less accurate than HL but 
better than HL1. We do not plot g%{r) to avoid obscuring 
the figures. 

VI. CONCLUSIONS 

Thus, the harmonic lattice model allows one to study 
static and dynamic properties of quantum and classical 
Coulomb crystals. The model is relatively simple, espe- 
cially in comparison with numerical methods like MC, 
PIMC and MD. The model can be considered as com- 
plementary to the traditional numerical methods. More- 
over, it can be used to explore dynamic properties of the 
Coulomb crystals and quantum effects in the cases where 
the use of numerical methods is especially complicated. 
For instance, the harmonic lattice model predicts singu- 
larities of the static inelastic structure factor at the posi- 
tions of Bragg diffraction peaks. We expect also that the 
HL model can describe accurately non-Coulomb crystals 
whose lattice vibration properties are well determined. 
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